import numpy as np
from EOS_collection import NRTL_MATRIX, antoine
import time

start_time = time.time()

# 常数
R = 8.314
atm = 101.325
Pa = atm
alpha = 0.3

# 基础参数 甲醇/水
a = np.array([1.201, -0.317])
b = np.array([1.001, 0.995])
Tb = np.array([[337, 373]])
x1 = np.array([[0.5]])
x2 = 1 - x1
x = np.concatenate((x1, x2))
T = np.sum(x.flatten() * Tb.flatten())
n = 12 #迭代循环次数 > 3

A = np.array([7.19736, 7.07404])
B = np.array([1574.99, 1657.46])
C = np.array([-34.29, -46.13])

dim = np.size(A)

def BUBBLE(Tb):

    for i in range(n):
        # 提供割线法的迭代点
        if i < 2:
            T = Tb[i]
        if i > 100:
            print(f"100次迭代不收敛！")
            exit(-1)

        # 计算活度系数
        tau_unknown = a + b / T
        tau = np.zeros((dim, dim))
        tau[np.where(np.eye(dim)==0)] = tau_unknown.flatten() #除了主对角线，其余依次赋值
        gamma = NRTL_MATRIX(x.T.reshape(-1, 1), tau, alpha).flatten() #将x变成列向量

        # 计算饱和蒸汽压
        Pvs = antoine(A, B, C, T)

        # 计算气相分率
        k = Pvs / Pa * gamma
        y = k.reshape(-1, 1) * x

        if i == 0:
            F0 = np.sum(y) - 1
            T0 = T
        else:
            F = np.sum(y) - 1
            # 满足迭代要求
            if np.abs(F) < 1e-5:
                break
            #割线法迭代
            T1 = T - F * ( T - T0 ) / ( F - F0 )
            T0 = T
            T = T1
            F0 = F
                
    return T


T = BUBBLE(Tb.flatten())
print(f"泡点温度是{T:.4f}K，{T-273.15:.4f}℃")
end_time = time.time()
elapsed_time = end_time - start_time
print(f"运行时间：{elapsed_time:.4f}秒")